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Extensive molecular dynamics simulations show that a short-range central potential, suited to 
model Ceo , undergoes a high temperature transition to a glassy phase characterized by the positional 
disorder of the constituent particles. Crystallization, melting and sublimation, which also take 
place during the simulation runs, are illustrated in detail. It turns out that vitrification and the 
mentioned phase transitions occur when the packing fraction of the system — defined in terms of an 
effective hard-core diameter — equals that of hard spheres at their own glass and melting transition, 
respectively. A close analogy also emerges between our findings and recent mode coupling theory 
calculations of structural arrest lines in a similar model of protein solutions. We argue that the 
conclusions of the present study might hold for a wide class of potentials currently employed to 
mimic interactions in complex fiuids (some of which of biological interest), suggesting how to achieve 
at least qualitative predictions of vitrification and crystallization in those systems. 



I. INTRODUCTION 

Central, short-range potential models of simple ana- 
lytic form have been the object of intense investigation 
over the last years. The main reason of such an inter- 
est is that they provide an approximate representation 
of the effective interactions and phase behavior experi- 
mentally observed in a variety of macrosized particle sys- 
tems, spanning from colloidal suspensions and protein or 
polymer solutions T, "5, 0, |^, la, B, Hi to other macro- 
molecular systems and fuUerenes O, 0, 0, 0| . 
A paradigmatic example in this respect is offered by pro- 
tein crystallization 0,0,^3- indeed, if one adopts a 
representation in which the macromolecules are spherical 
particles interacting through a central short-range po- 
tential, and the solvent is assumed to be a structureless 
continuum, the liquid-vapor binodal and the sublimation 
line thereby calculated exhibit all features of the protein- 
rich/protein-poor and crystallization line observed in real 
globular protein solutions. These models then offers an 
unvaluable tool for the study, by usual means of statisti- 
cal mechanics and simulation, of the relative location of 
phase coexistence lines in real systems. 

This scenario has been further enriched by recent in- 
vestigations of glass transition in one of the most stud- 
ied short-range model, namely the attractive hard-core 
Yukawa fluid. In this context the mode coupling theory 
predicts the onset of two glassy phases deeply inside the 
metastable liquid-solid region of the model LZ|. The for- 
mation of distinct vitreous phases with peculiar internal 
characteristics has also been reported in recent experi- 
ments on colloid-polymer mixtures [TsL IT^ l20j . 

These evidences prompted us to an extensive investiga- 
tion of phase coexistence and glass transition conditions 
in short-range models, with particular attention to the 
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rich phenomenology which characterizes the metastable 
region enveloped between the freezing and the melting 
line of such systems. The paper focuses in particular 
on the Girifalco potential for fuUerene Cgo f2l|. This 
model can be considered as "marginally" short-range in 
the sense that, whereas in colloidal and protein systems 
the interaction between (model) macroparticles reduces 
practically to zero within a fraction of the particle diam- 
eter and the liquid-vapor equilibrium is only metastable, 
in the Girifalco model the decay length of interactions is 
slightly greater than the fullerene diameter and a stable 
liquid phase does exist, albeit confined to a temperature 
interval smaller than 100 K ^9.. .12„.13».14) |. with no experi- 
mental evidence of it hitherto reportedjlllllillllllil . 
Such characteristics appear consistent with early spec- 
ulations about the "elusive diffusive" nature of liquid 
Ceo jOl > as well as with recent evidences of the absence, 
in this peculiar fluid, of cage effects in the velocity auto- 
correlation functions, and of collective phenomena, ordi- 
narily observed in simple liquids p7| . 

We have recently reported in a preliminary communi- 
cation (hereafter referred to as I) Molecular Dynam- 
ics (MD) evidences of the onset of a glass transition in 
Cgo at T = HOOK upon a pressure of 3.5 MPa, as as- 
sociated to positional disorder. Such a glassy phase is 
quite distinct from the experimentally detected Cgo ori- 
entational glass, which forms at 90 K as a consequence of 
the freezing at low temperature of the residual disorder 
in the orientation of the fullerene cages • 

As preliminary documented in I, vitrification of Cgo at 
high temperature is assessed through quenching cycles of 
the liquid phase at various pressures, by analyzing a num- 
ber of thermodynamic, structural and dynamical quanti- 
ties. Here, we report and discuss in detail the MD results 
we cumulated to support our conclusion that a positional 
glass of Ceo is effectively formed. In the same context, 
we discuss other transitions undergone by the system in 
the course of cooling or heating cycles, such as crystal- 
lization, melting and sublimation. As we shall detail. 
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such transitions take place deeply inside the metastable 
region in coincidence of the crossing of two density-vs- 
temperature loci over which the "effective" packing frac- 
tion of the system equals that of hard spheres at their 
melting and glass transition, respectively. The analysis 
of such hard-sphere-like behavior enlights the analogies 
of our study with the case, discussed in reference of 
the onset of glassy phases in model protein solutions, sug- 
gesting that our results can be useful in the more general 
context of complex fluids investigation. 

The paper is structured as follows: in section II we in- 
troduce the model and describe the computational strate- 
gies. In III results for thermodynamics, structural and 
dynamical properties are reported. The discussion and 
conclusions follow in section IV. 



eventually obtained at the end of the cooling cycles, or 
from a perfect fee crystal at room temperature with 864 
particles. Other details are given in the presentation of 
results. 



III. RESULTS 

For clarity sake, we shall detail in the first instance the 
transitions undergone by the system during slow cooling 
or heating sequences, followed by a presentation of results 
obtained during the quenching routines. Just before, we 
introduce the concept of "effective" packing fraction for 
the system at issue, with two associated loci of thermo- 
dynamic states. 



II. MODEL AND SIMULATION STRATEGIES 



We study a system composed of particles interacting 
via the Girifalco potential |2]| . 
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where s = r/d, ai = N^A/12d^, and a2 = iV^S/QOd^^. 

= 60 and d = 0.71 nm are the number of carbon atoms 
and the diameter, respectively, of the fuUerene particles, 
A = 32 X 10"^*' ergcm^ and B = 55.77 x IQ-^^^ ergcm^^ 
are constants entering the 12-6 potential (f>{i')— ~A/r^ + 
B /r^^ through which two carbon sites on different spher- 
ical molecules are assumed to interact. The distance 
where the potential crosses zero, the position of the 
potential well minimum and its depth, are a = 0.959 nm, 
7'min = 1.005 nm, and e — 0.444 x 10~i2erg, respec- 
tively |2l|. We assume i;(r) = oo for s < 1. 

We have used the MD method introduced by Ander- 
sen [sj to simulate constant-pressure, constant-enthalpy 
(NPH) conditions, and the Verlet algorithm to integrate 
the equations of motion over time steps Ai = 5 fs; cu- 
bic boxes and periodic boundary conditions are assumed. 
Most results are obtained for N = 864 and N = 1000 to- 
tal number of particles; the smaller sample is fully com- 
patible with an initial fee lattice arrangement of the par- 
ticles in a cubic simulation box, whereas the 1000-particle 
distribution obviously implies defects. Both sizes are part 
of a more general analysis employing also five hundred 
and 1372 particles, to exclude any dependence of simu- 
lation results on the sample size. A wide range of pres- 
sures, spanning from close-to-atmospheric conditions up 
to 250 MPa (2.5kbar), has been explored to verify the 
effects on crystallization and glass transition conditions. 
Cooling and quenching cycles start from an initial super- 
critical or liquid configuration (see figure^, identified as 
such on the basis of the known phase diagram of Ceo llJl ■ 
Heating cycles begin either from the defective crystals 



A. Effective hard sphere diameters and 
constant "packing" loci 

We introduce an effective hard sphere diameter for 
the Cgo-Cgo interaction with the definition of a refer- 
ence potential for v(r) of equation (^) according to the 
well-known Weeks-Chandler-Andersen (WCA) prescrip- 
tion 0, 
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and adopt the Barker and Henderson expression for the 
effective hard-core diameter 331, 
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The Barker-Henderson prescription for the hard sphere 
diameter is by no means unique. Other more sophisti- 
cated recipes, hinging for instance on the use of the cavity 
distribution function y{r) = exp[/3'i;(r)](7(r) [s^l (where 
g(r) is the radial distribution function), have been pro- 
posed in the literature. We have verified however that 
the resulting changes in the estimate of the effective di- 
ameter with respect to the simpler BH prescription are 
quite minor. 

We now look for densities p = p(r) which make the 
effective packing fraction r\ = ir/Gpcr'^YL constant with 
respect to temperature variations, and equal to a pre- 
fixed value. In particular, we impose rj to be equal to 
the packing of hard spheres at their own melting and 
glass transition, respectively. The value 77^^^ = 0.545 has 
been reported for the packing fraction of hard spheres at 
melting [s^. A packing of 0.55 has also been reported 
in constant pressure MD simulation of hard spheres . 
Here we choose an intermediate value (see I), by imposing 



7r/6pm(T)cr|H 



0.548. 



(4) 
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As for the packing of hard spheres at the glass tran- 
sition, 1]^^ , the value 0.58 has been reported in refer- 
ence [23 ■ whereas effective packing fractions at vitrifica- 
tion for continous potentials, mostly of Lennard-Jones- 
likc form (see, e.g. |^ and references), approach 0.58 
with a considerable dispersion of values due to the differ- 
ent prescriptions either for the effective diameter or for 
the determination of the glass transition temperature. 
We require in particular (see I): 

r?g = 7r/6pg(r)aiH = 0.574. (5) 

On the basis of equations l@J and Q we can determine 
the two functions Pm{T) and Pg(T) which are shown in 
figure n As we shall see, these two p{T) loci correlate 
significantly with the crystallization and vitrification con- 
ditions of the Ceo model. 

B. Cooling cycles and crystallization; heating 
cycles and sublimation or melting 

We first illustrate the case when the pressure P = 
3.5 MPa. Cooling starts from a high temperature fluid 
configuration at T = 2000 K (see figure ^ top). In se- 
quential stepwise drops AT = 30 K, the system evolves 
over 20 000 steps (correspondiong to 100 ps) at fixed tem- 
perature, followed by 10 000 cumulation steps where the 
system evolves freely. The final temperature T, box vol- 
ume V, enthalpy H, and other thermodynamic, struc- 
tural and dynamical properties are recorded at each AT 
(see below) . Cooling is always arrested when room tem- 
perature is achieved or sligthly before this threshold. Sta- 
tistical uncertainties on T, p, and P turn out to be 0.8 K, 
0.0003 nni^'^ and 0.1 MPa, respectively. 

Cooling paths are visualized in figuren(top) in the p-T 
representation of the phase diagram of Cgo . Volumes and 
enthalpies changes with the temperature are displayed in 
figure 121 (top) ; as visible (and already shown in I for the 
case N = 1000), for both system sizes investigated V 
and H undergo a marked drop at T ~ 1300 K, accompa- 
nied by a temporary increase of the temperature. Such 
a highly nonmonotonic behaviour of the two thermody- 
namic quantities is fully consistent with a transition of 
the system to the solid phase as further documented by 
(i) the behavior of the radial distribution function, g(r) 
for T < 1307 K (figure 13 top), (ii) the temperature de- 
pendence of the diffusion coefficient D (figure ^top) and 
(Hi) the snapshots of the Ceo configurations at differ- 
ent temperatures (figure IS)). The incipient solid phase is 
characterized by a fee crystalline arrangement, as visible 
from an inspection of the peaks' positions in the g{r). 

The variation of N from 1000 to 864 docs not produce 
any significant effect either on the overall cooling path, or 
on the location of the turning points which marks the on- 
set of crystallization. Such an outcome is obtained if we 
evolve the smaller system for a longer time. In fact, under 
equal time elapsed conditions the turning point with 864 
particles would fall at lower temperatures. We interpret 



TABLE L Crystallization (c), melting (m) and glass (g) tran- 
sition parameters. Melting data refer to the defective crystal 
with TV = 1000. Pressures are given in MPa, temperatures in 
K, densities in nm~"^, and diameters in nm; in the last column 
the packing fraction rj^ — 7r/6pxO"|jj. 

P Tc pc ffBH Vc 

3.5 1307 1.124 0.976 0.547 

40 1464 1.130 0.975 0.548 

150 1880 1.138 0.972 0.547 

250 2200 1.150 0.970 0.550 

3.5 1931 1.144 0.9716 0.549 

P Tg Pg Q-BH Ve 
3.5 1100 1.168 0.978 0.572 

40 1170 1.177 0.977 0.575 
150 1480 1.187 0.975 0.576 
250 1700 1.190 0.973 0.574 



this result as a manifestation of the fact, also emerging in 
heating cycles, that the 1000-particle case, characterized 
by a larger simulation box and defects with respect to 
the perfect fee arrangement, easier allows for the inter- 
nal particle rearrangements which trigger the phase tran- 
sitions. We also note that the first of the turning points, 
signalling the onset of crystallization falls practically in 
coincidence of the crossing between the cooling path and 
the 1]^^ = 0.548 locus introduced in equation Q (com- 
pare figure n^^nd figure EJ; we shall discuss this point in 
the next section. 

No appreciable effect on the system behavior is exerted 
by the cooling rate since cooling cycles with AT = 15K 
lead to substantially similar results (see I). Sensitive 
changes are instead associated to the value of the im- 
posed pressure: the modifications undergone by the cool- 
ing paths as P varies from 3.5 to 150 and 250 MPa are 
shown in figure ^ (middle panel, see also the inset) and 
figure 121 (compare top and middle panels); the radial dis- 
tribution functions for the case P = 250 MPa are shown 
in the bottom panel of figure and numerical values of 
the crystallization parameters are reported in tabled It 
appears that the increase of the pressure rises dramati- 
cally the temperature where the onset of crystallization 
occurs. The densities attained at crystallization, and at 
the end of the cooling paths, are close to each other and 
exhibit a trend to increase with the pressure. These den- 
sities are all smaller than those of the perfect crystal 
heated at the corresponding temperature, indicating that 
defective crystals are formed. One can realize that this is 
so by comparing the physical characteristics of such crys- 
tals with those typical of a perfect crystal of 864 particles. 
We consider to this aim crystals at room temperature as 
previously obtained through the cooling from high tem- 
peratures of samples with N — 1000 and N = 864, and 
a perfect crystal of 864 particles under the same temper- 
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ature and pressure conditions. As visible from figure ^ 
and figure 12 the final densities of the samples obtained 
from the cooling procedure are smaller than that of the 
perfect 864 crystal. Such lower density values must be 
attributed to the presence of voids in the crystalline ma- 
trix formed during the cooling process. The comparison 
of radial distribution functions in figure |S1 enlights the 
smoother structure of cooled samples with respect to the 
perfect crystal. At P=250MPa, after the onset of crys- 
tallization, the cooling path runs almost overlapped with 
the melting curve of Cgo Q (see figure Q middle). The 
high pressure is thus able to force the system to evolve 
out of the metastable region approximately along the true 
coexistence line. 

The transition to the crystal takes place under strong 
supercooling conditions, well beyond the freezing line, 
indeed, when the latter is crossed, no feature in the radial 
distribution function signals the onset of any structural 
order typical of the solid phase (see figure OJ. Signatures 
albeit faint of such a crossing can be recognized however 
in the V and H patterns which show a corresponding tiny 
nonmonotonicity (see top panels of figure |5| and insets), 
as well as in the diffusion coefficient (see figure 01 and 
insets). As visible, runs with different number of particles 
exhibit the same features, ruling out the possibility that 
our observations are merely a consequence of statistical 
noise. 

A similar situation emerges when the systems ap- 
proaches the metastable portion of the liquid-vapor bin- 
odal of Cgo, calculated in reference [s^ and reported in 
figure n It appears that the cooling paths get closer 
to the binodal the lower the pressure becomes. At 
P ~ 3.5 MPa a crossing can be extrapolated to take 
place at T = 1700 K and p — 0.95 nm~'^, a state point 
where V and H patterns show another nonmonotonicity, 
magnified in the insets of figure |2 Remarkably, also the 
diffusion coefficient in figure 0] is nonmotonic in corre- 
spondence of the estimated crossing of the binodal. We 
have carried out very long simulation runs (up to sev- 
eral million time steps) in correspondence of the crossing 
of the freezing and binodal lines, to check whether the 
system is able to develop any signal of the incipient tran- 
sitions. In no case, however, we could monitor transfor- 
mations in the thermodynamic and diffusion coefficient 
behavior comparable to those heralding at the crystalliza- 
tion temperature T = 1307 K, as previously described. 
Several similarities emerge in comparison with the sim- 
ulation study that two of us have performed |4Clj| on a 
modified Lcnnard-Jones potential used to model globu- 
lar protein solutions 0]. In that case, the transition to 
the crystalline arrangement takes place during isochoric 
coolings of the system deeply beneath the freezing line, 
when the system reaches the metastable binodal line. In 
the present case, the system must be cooled even substan- 
tially below the metastable binodal before crystallization 
can start. 

Starting from the two solid simulation samples eventu- 
ally obtained through the cooling cycles with N = 1000 



and 864, and from a third one built as a perfect crystal 
with 864 particles, we gradually heat all these systems 
through successive AT=30K increases, up to the tem- 
perature at which they undergo an abrupt transition to 
a fluid configuration. Bottom panels of figures and El 
refer to this sequences; it appears that the defective sam- 
ple with one thousand particles undergoes a true melt- 
ing since the systems jumps to a thermodynamic point 
inside the liquid pocket of the phase diagram whereas 
the 864 particle crystal, either defective or perfect, is 
instead substantially overheated until it jumps directly 
to a vapor phase, thus exhibiting sublimation. The role 
of defects at intermediate temperatures emerges by the 
comparison of the radial distribution functions obtained 
from the heating of the perfect crystal at relatively high 
temperature with the corresponding g{r) of a crystal at 
an equal temperature obtained from cooling: in figure 
one can see that all features are less sharp in the latter 
case, with some peaks even missing with respect to the 
perfect crystal case. As can be seen in the corresponding 
panels in figure (21 an evident "hysteresis" characterizes 
the heating cycles with respect to the cooling sequences. 



C. Quenching and glass transition 

Quenching of the liquid at different pressures is car- 
ried out for N = 1000 and 864, starting from an ini- 
tial temperature T — 1950 K, either through a sequence 
of AT = 150 K decrease steps, or through a sudden 
AT=1000K drop. Temperature variations are imposed 
over 20 000 simulation steps, followed by 10 000 steps of 
free evolution. The quenching paths are displayed in fig- 
ure n (top and middle panels) while the volume and en- 
thalpy behavior is visible in figure (21 (top and middle 
panels). 

The onset of a glassy phase during both quenching pro- 
cedures is visually documented in the snapshots shown in 
figure [31 and is proved by several evidences on thermody- 
namic, structural and dynamical quantities. We first ob- 
serve that the V and H patterns in figure [3 show a glass 
branch running above the crystallization line; moreover, 
at variance with what we have observed during the slow 
cooling, the radial distibution functions in figure do 
not develop the peak structure of the fee arrangement 
at low temperatures; rather, down to ~ 300 K they ex- 
hibit the twofold structure of the second peak typical 
of the glass. Such a structural evidence has a counter- 
part in the diffusion coefficient behavior which, as visi- 
ble in figure 01 does not exhibit the drop associated to 
the onset of crystallization process, but manifests only a 
change of slope. On the other hand, the mean square dis- 
placements show a plateau for temperatures T < HOOK, 
signalling structural arrest, also confirmed by the shape 
variations of the velocity autocorrelation function with 
the temperature (see figure [TJ. Changes of slope as the 
temperature decrease, typically associated to the glass 
transition, occur also in V and H (figure O, as well as 



5 



in the thermal expansivity a, the specific heat Cp, and 
the Wendt- Abraham ratio R of the first peak to the first 
minimum height in g{r) (figure [HJ- 

Evidence of the glass formation also comes from the 
shear viscosity which we have calculated via the Green- 
Kubo relation lifl: 
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is the off-diagonal component of the stress tensor, vf and 
Xij represent, respectively, the x component of the ve- 
locity and separation distance, r.^, between molecules' 
center-of-mass i and j, and fy is the y component of the 
force exerted on atom i by atom j. As can be appreci- 
ated in figure El 'ysh exhibits a marked increase across the 
glass transition region and displays an Arrhenius behav- 
ior as a function of T (see inset), indicating that the glass 
formed is a "strong" one 42]; the orientational glass of 
Geo is similarly "strong" [29t..3Q]. 

The glass transition parameters are obtained from the 
intersection of the two (extrapolated) branches in the 
V, H, Cp, a and R patterns as functions of the tem- 
perature (see reference ^3 details): we estimate at 
P = 3.5 MPa Tg ~ HOOK and Pg ~ 1.168 nm'^. 

As visible from figure |2 a glass transition is observed 
also at pressures higher than 3.5 MPa, at the cost to in- 
crease the quenching rate with respect to 1.5 x lO^'^ K/s 
adopted at 3.5 MPa. For instance, at P = 250 MPa no 
glassy phase is obtained if AT < 300 K over 20 000 time 
steps. Glass transition temperatures and densities at dif- 
ferent pressures are reported in table It is immediate 
to verify that, similarly to what observed with the crys- 
tallization transitions, the glass transitions (at different 
pressures) occur for states falling at the crossing of the 
quenching paths with the p{T) locus defined by equa- 
tion lO, where 77g= ry^^. We shall comment on these 
evidences in the next section. 

The quenching rates we adopt in our simulations are 
much higher than those presently achieved in experi- 
ments (typically 10^ K/s). We have however prelimi- 
nary simulative evidences that in equimolar Girifalco 
Ceo/Cro and Cgo/Ggg mixtures crystallization does not 
occur even upon cooling rate one or two order of mag- 
nitude lower than those adopted here; the Ggo/Gyo mix- 
ture, for instance, remains liquid during the cooling pro- 
cedure down to HOOK, a temperature lower than that 
at which amorphization of pure Cgo fuUerite occurs p3 | 
(see also note 26j). These findings make it plausible that 
a glassy phase might be formed experimentally at least 
from mixed fuUerene systems. 



IV. DISCUSSION AND CONCLUSIONS 

The effective packing at crystallization for different 
pressures have values close to the packing fraction of hard 
spheres at melting ?7j„^=0.545, as reported in tabled As 
a matter of fact, it thus apperas that the rjax{T) locus 
introduced in equation Q acts as an interpolation curve 
among the transition points determined a posteriori, and 
can be used to forecast crystallization temperatures and 
densities at pressures other than those here considered. 
A similar evidence is found for the effective packing at 
the glass transition ?7g(T) of equation jSJ also reported 
in the table, that is quite close to the hard sphere glass 
transition value, 77^^ = 0.58 [s^]. As for the onset of crys- 
tallization and glass transitions, an obvious deduction is 
that the Girifalco model follows rather closely a hard- 
sphere-like behavior. Our evidences, however, allow us to 
get a deeper insight into the role of the Pm{T) and Pg{T) 
loci. Indeed, we have found that when supercooling is 
pushed beyond Pm{T) the system invariably crystallizes; 
by contrast, a state close to, but on the left of pa\{T), 
remains in the supercooled phase even over the longest 
simulation run we could perform, namely 12 million sim- 
ulation steps, corresponding to 60 ns. It thus appears 
that the Pm{T) locus approximately corresponds to an 
instability boundary of the metastable region. This ob- 
servation also emerges from Monte Carlo calculations of 
the free energy of the Girifalco model at T = 2100 K 
(see inset of figure^; in fact, pm(2100K) coincides with 
the point where the free energy shows an inflection point 
associated with a sudden drop of the pressure of the sim- 
ulation sample, i.e. with a mechanical instability of the 
system. A similar property might characterize the free 
energy branch generated by approaching the metastable 
region from the solid side and since the interval sepa- 
rating the two branches in figure ^ is quite narrow, it 
is conceivable that our simulation strategy is unable to 
discriminate between the densities at the crystallization 
and melting instabilities. It is interesting to observe that 
in Molecular Dynamics simulations of the Lennard- Jones 
potential, crystallization occurs much beyond the melting 
line, deeply inside the solid region Indeed we have 
carried out MD simulations of the same potential and 
verified that states located inside the metastable fluid- 
solid region, do not undergo crystallization even after 35 
million time steps (~ 180 ns). 

Our results fot the glass transition positively agree 
with the predictions made by Foffi and coworkers Q on 
the glass transition in protein solutions modeled through 
the attractive hard-core Yukawa (HGY) potential. These 
authors find that for inverse decay length z of the Yukaw- 
ian term spanning from 5 to 60 times the reciprocal of the 
particle diameter, the glass transition line is an almost 
vertical locus with a T = 00 asymptote i] = i]^^ = 0.58. 
On the other hand, as discussed in [43 , the physical prop- 
erties of the Girifalco model can be reasonably repro- 
duced by a Yukawa potential with z ~ 4 (though with 
moderate differences in the phase diagram), close indeed 
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to the lowest value investigated in reference [7| . A quali- 
tative agreement thus emerges between the prediction of 
the glass transition in the two different systems. In the 
HCY model the formation of "repulsive" and "attractive" 
glasses depends on the value of z |7j ; since at z as low as 
five only the repulsive glass would be formed, we identify 
our Ceo "positional" glass as a repulsive one. 

In conclusion, we have documented that effective hard- 
core exclusion effects play an important role for the deter- 
mination of the glass transition in a short-range potential 
well suited to model various fuUerenes. Mode Coupling 
Theory calculations for the onset of the glass transition 
in a short-range model of globular protein solutions lead 
to similar conclusions at least for relatively low values of 
the decay potential parameters. In this context, it could 
be worth to carry out specific Mode Coupling calcula- 
tions for the present model fuUerene in order to compare 
the structural arrest line thereby predicted with our re- 



sults for the glass line based on simulations and Pg(T) 
locus predictions. Moreover, colloidal suspensions and 
colloid-polymer mixtures show a behavior characterized 
by several similarities with the observations reported in 
this work 6] . We argue that the basic equation Q might 
hold to a high accuracy also for other interaction poten- 
tials, and hence it could provide a simple framework for 
a qualitative but immediate prediction of the glass line 
in a variety of model systems. 
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FIG. 1: Top: Continuous lines: phase diagram of Ceo [l^ with superimposed the Pm(r) [long-dashed line, equation ((IJ] and 
Pg(T) [short-dashed line, equation ||KJ] loci introduced in this work. Circles and squares: cooling with 1000 and 864 particles 
at P = 3.5 MPa, respectively. Inset: Helmholtz free energy vs p at T = 2100 K [ijj- Middle: Effect of increasing pressure on 
the cooling and quenching paths. Full circles, squares and triangles: cooling paths with A'^ = 1000 at pressure P = 3.5, 150 and 
250 MPa, respectively; open circles, squares and triangles: quenching with A'^ — 1000 at P = 3.5, 150 and 250 MPa ; di amonds: 
quenching at 3.5 MPa with 864 particles. The dot-dashed line indicates the metastable portion of the binodal [39|. Inset: 
crystallization density as a function of the pressure and corresponding temperatures. Bottom: Heating paths of the defective 
and perfect crystals. Circles and squares: heating of the defective crystal with A'' — 1000 and 864, respectively; diamonds and 
triangles: heating of the perfect crystal with — 864 at P = 3.5 MPa and 12 MPa, respectively. 
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FIG. 2: Top: Volume (left) and enthalpy (right) as functions of the temperature in cooling and queneliiug cycles at P = 3.5 MPa. 
Full circles, open triangles: cooling with N = 1000 and 864, respectively; open circles, Eisterisks: quenching with A'^ = 1000 and 
864, respectively. Inset: expansion of the crossing zones (encircled) of the freezing and binodal line during cooling. Middle: V 
(left) and H (right) vs T in cooling (full circles) and quenching (open circles) cycles with A'" = 1000 at P = 250 MPA. Bottom: 
V (left) and H (right) vs T in heating cycles. Open circles, full diamonds: heating of the defective crystal with A' = 1000 
and 864 at P = 3.5 MPa, respectively; open diamonds, full triangles: heating of the perfect crystal with A^ = 864 particles at 
P = 3.5 and 12 MPa,respectively. Inset: expansion of the volume scale around the transition to the liquid (for A^" = 1000 at 
3.5 MPa, open circles), and the sublimation to a highly expanded vapor phase (for all other cases). 
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FIG. 3: Top: evolution of g{r) with the temperature during cooling at 3.5 MPa, with N = 1000 (full line) and N = 864 
(dashed lino). Temperatures aside the curves refer to the — 1000 (left) and 864 (right). Middle: same as in the top panel 
for quenching cycles. Bottom: Cooling (full line) and quenching (dashed line) for N = 1000 at P = 250 MPa. 
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FIG. 4: Top: Diffusion coefficient of cooling (open circles, full triangles) and quenching (open diamonds and full squares) 
cycles for N = 1000 and 864, at 3.5 MPa, respectively. Upper inset: linear scale expansion of the D behaviour at the crossing 
of the freezing and binodal lines; lower inset: Arrhenius plot of the diffusion coefficient in quenching cycles. Bottom: Diffusion 
coefficient in cooling (full triangles) and quenching (open circles) cycles at P = 250 MPa for N = 1000. Inset: Arrhenius plot 
of D in the quenching cycles. 
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FIG. 5: Snapshots of the A'' = 1000 sample during the coohng (top) and the quenching (bottom) sequence at P = 3.5 MPa 
and T = 1937, 1307, 1152, and 300 K (from left to right). 
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FIG. 6: Evolution of g{r) during the heating sequence for the perfect crystal with N = 864 (full line) at T = 598, 1030, 
2346, and 2391 K (from top to bottom), and for the defective crystal with iV = 1000 (dotted line) at T = 584, 1037, 1931, and 
1934 K. The dashed line is g{r) at T — 1022 K, as obtained by cooling the A'^ — 864 sample from high temperature (1937 K). 
The comparison with the 1037 K heating case enlights the smoother features resulting from cooling. 




FIG. 7: Top: Mean square displacement (in A/s units) in the A'^ = 864 sample as a function of the simulation time (in 
picoseconds) during the quenching sequence at P = 3.5 MPa. Bottom: Evolution of the velocity autocorrelation function 
during the same sequence (to=time step). 




FIG. 9: Shear viscosity vs temperature during a quenching sequence of the N = 864 sample at P = 3.5 MPa. The inset shows 
the Arrhenius behaviour in the approach to the glass temperature. 



